Incremento del nivel medio del mar

Juan Sebastián Mendoza Páez

Estudiante Ing. Forestal

¿Por qué es importante?

NASA

“Desde 1880, el nivel del mar global ha aumentado 20 centimetros. Para el 2100 se proyecta que aumente entre 30 y 122 centimetros más.”

Bases de datos

Fuentes

  • Oficina Nacional de Administración Oceánica y Atmosférica (NOAA)
  • Organización de Investigación Científica e Industrial del Commonwealth (CSIRO)
Mostrar Código 
pacman::p_load(TSstudio, tsbox, tsoutliers, tidyverse, equatiomatic,
               lubridate, xts, magrittr, forecast, dygraphs)

url <- 'https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt'
cov <- read.table(file = url) %>%
  select(3, 4) %>%
  rename(fecha = 1, co2 = 2) %>%
  mutate(fecha = date_decimal(fecha) %>%
           floor_date(unit = 'month')) %>%
  filter(fecha <= '2013-12-01')

data <- read.csv('gmsl.csv') %>%
  mutate(fecha = date_decimal(Time) %>%
           floor_date(unit = 'month')) %>%
  rename(gsml = 2) %>%
  select(4, 2)%>%
  filter(fecha > '1958-02-01') %>%
  left_join(y = cov) %>%
  {xts(x = .[,2:3], order.by = .$fecha)}

rm(cov)

dygraph(data, main = 'Nivel medio del mar global vs
        Concentración de CO2 en la atmósfera', height = 400) %>%
  dySeries("co2", axis = 'y2', label = 'CO2') %>%
  dySeries('gsml', label = 'GSML') %>%
  dyAxis('y', label = 'Nivel del mar (mm)') %>%
  dyAxis('y2', label = 'CO2 (ppm)', independentTicks = T) %>%
  dyRangeSelector()

Estacionalidad y prueba de raíces unitarias

Mostrar Código 
acf(train[,'gsml'], main = '')

Mostrar Código 
tseries::adf.test(x = train[,'gsml'])

    Augmented Dickey-Fuller Test

data:  train[, "gsml"]
Dickey-Fuller = -2.7778, Lag order = 8, p-value = 0.2491
alternative hypothesis: stationary

Modelos

Mostrar Código 
matriz_diseño <- model.matrix(~ -1 + train[,'co2'])
colnames(matriz_diseño) <- 'co2'

modelo <- Arima(y = train[,'gsml'], 
                order = c(2,1,3),
                xreg = matriz_diseño,
                include.drift = T)

modelo %>% lmtest::coeftest()

z test of coefficients:

       Estimate Std. Error  z value  Pr(>|z|)    
ar1    0.797070   0.054991  14.4947 < 2.2e-16 ***
ar2   -0.095533   0.054919  -1.7395   0.08194 .  
ma1   -0.339090   0.040237  -8.4274 < 2.2e-16 ***
ma2    0.208116   0.038717   5.3753 7.646e-08 ***
ma3   -0.693544   0.028246 -24.5535 < 2.2e-16 ***
drift  0.208518   0.041471   5.0281 4.955e-07 ***
co2   -0.171228   0.087716  -1.9521   0.05093 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC: 2595.726
Mostrar Código 
mod <- Arima(y = train[,'gsml'], 
                order = c(2,1,3),
                include.drift = T)

mod %>% lmtest::coeftest()

z test of coefficients:

       Estimate Std. Error  z value  Pr(>|z|)    
ar1    0.794860   0.055187  14.4031 < 2.2e-16 ***
ar2   -0.092967   0.055205  -1.6840   0.09218 .  
ma1   -0.341537   0.040571  -8.4182 < 2.2e-16 ***
ma2    0.204242   0.039030   5.2329 1.668e-07 ***
ma3   -0.692204   0.028275 -24.4809 < 2.2e-16 ***
drift  0.187589   0.039109   4.7965 1.614e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC: 2597.526
Mostrar Código 
holt <- holt(y = train[,'gsml'], h = 12, initial = 'optimal')
holt$model
Holt's method 

Call:
 holt(y = train[, "gsml"], h = 12, initial = "optimal") 

  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.6589 

  Initial states:
    l = -49.0075 
    b = -2.4967 

  sigma:  2.4035

     AIC     AICc      BIC 
5429.938 5430.030 5452.384

Modelos

Mostrar Código 
extract_eq(modelo, wrap = T, use_coefs = T)

\[ \begin{alignat}{2} &\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +0.21\operatorname{t} -0.17\operatorname{co2}_{\operatorname{t}} +\eta_{t} \\ &\text{where}\quad &&(1 -0.8\operatorname{B} +0.1\operatorname{B}^{\operatorname{2}} )\ (1 - \operatorname{B}) \eta_{t} \\ & &&= (1 -0.34\operatorname{B} +0.21\operatorname{B}^{\operatorname{2}} -0.69\operatorname{B}^{\operatorname{3}} )\ \varepsilon_{t} \\ &\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})} \end{alignat} \]

\[( 1 \ - \ 0.79B \ + \ 0.09B^2) \ \cdot \ (1 - B) \ \cdot \ (y_t - 0.19t) \\ = (1 \ - \ 0.34B \ + \ 0.2B^2 \ - \ 0.69B^3) \ \cdot \ \epsilon_t\]

Supuestos

Mostrar Código 
mod %>%
  checkresiduals(theme = theme_bw(), test = F)

Mostrar Código 
mod %>%
  residuals() %>%
  shapiro.test()

    Shapiro-Wilk normality test

data:  .
W = 0.99855, p-value = 0.8733
Mostrar Código 
mod %>%
  residuals() %>%
  Box.test(type = 'Ljung-Box')

    Box-Ljung test

data:  .
X-squared = 0.028529, df = 1, p-value = 0.8659

Outliers

Mostrar Código 
tso(y = train[,'gsml'], xreg = train[,'co2'], tsmethod = 'auto.arima', args.tsmethod = list(stepwise = F, approximation = F))
Series:  
Regression with ARIMA(2,1,3) errors 

Coefficients:
         ar1      ar2      ma1     ma2      ma3   drift     xreg
      0.7971  -0.0955  -0.3391  0.2081  -0.6935  0.2085  -0.1712
s.e.  0.0550   0.0549   0.0402  0.0387   0.0282  0.0415   0.0877

sigma^2 = 2.993:  log likelihood = -1289.86
AIC=2595.73   AICc=2595.95   BIC=2631.63

No outliers were detected.

Predicciones

Mostrar Código 
prediccion <- forecast(object = modelo, h = 12, level = 95, xreg = test[,'co2'])

serie <- prediccion %>%
  as_tibble() %>%
  ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>%
  {cbind(data[, 'gsml'], .)}

colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr')

serie %>%
  dygraph(height = 350) %>%
  dySeries('gsml', label = 'Observada') %>% 
  dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>%
  dyRangeSelector()
Mostrar Código 
accuracy(prediccion, test[,'gsml'])[,c(2,5)]
                 RMSE     MAPE
Training set 1.719463      Inf
Test set     9.180928 12.04332
Mostrar Código 
serie <- holt %>%
  as_tibble() %>%
  select(1, 4, 5) %>%
  rename(predico = 1, lwr = 2, upr = 3) %>%
  ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>%
  {cbind(data[, 'gsml'], .)}

colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr')

serie %>%
  dygraph(height = 350) %>%
  dySeries('gsml', label = 'Observada') %>% 
  dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>%
  dyRangeSelector()
Mostrar Código 
accuracy(holt, test[,'gsml'])[,c(2,5)]
                 RMSE     MAPE
Training set 2.396217      Inf
Test set     4.513079 5.701727

Gracias

Series de Tiempo Univaridas